
# Packages ----------------------------------------------------------------

library(tidyverse)
library(haven)
library(sf)
library(rmapshaper)
library(tmap)


# Setup and load data -----------------------------------------------------

# Path to project directory
proj_dir <- Sys.getenv(c("Rep_smokelabor"))

# Coordinate reference systems

# - NAD83 geodetic coordinate system: see https://epsg.io/4269
# Used by the United States Census TIGER
NAD83 <- "+init=epsg:4269"

# - Albers projected coordinate system: see https://epsg.io/3083
# Used by the the United States Census Bureau
Albers_US <- "+init=epsg:3083"

# US shape files
US_states <- file.path(proj_dir, "1_build/geo/data/tiger_simple/state/2010/ts_2010_us_state_05pct.rds") %>%
  readRDS() %>%
  st_transform(Albers_US) %>% 
  ms_simplify(keep = 0.25, keep_shapes = TRUE)

US_counties <- file.path(proj_dir, "1_build/geo/data/tiger_simple/state/2010/ts_2010_us_county_05pct.rds") %>% 
  readRDS() %>% 
  st_transform(Albers_US) 
  # ms_simplify(keep = 0.04, keep_shapes = TRUE) 

# County data to overlay on the US map
df <- "1_build/hms/raw/county_smoke_annual.dta" %>% 
  file.path(proj_dir, .) %>% 
  read_dta() %>% 
  rename(GEOID10 = COUNTY10)

US_counties <- US_counties %>% 
  mutate(state_group = case_when(
    STATEFP10 == "02" ~ "Alaska",
    STATEFP10 == "72" ~ "Puerto Rico",
    STATEFP10 == "15" ~ "Hawaii",
    TRUE ~ "Continental US")) %>% 
  rename(state = STATEFP10, 
         county = COUNTYFP10)

# Join with estimations, limit to contiguous US states
US <- inner_join(US_counties, df, by = c("GEOID10")) %>% 
  filter(!state_group %in% c("Alaska", "Hawaii", "Puerto Rico"))

# Interior state borders for contiguous US
US_int_borders <- US_states %>% 
  filter(!NAME10 %in% c("Alaska", "Hawaii", "Puerto Rico")) %>% 
  st_geometry() %>% 
  st_cast(to = "MULTILINESTRING") %>% 
  st_difference(st_cast(st_union(US_states), to = "MULTILINESTRING")) %>% 
  st_cast("MULTILINESTRING")

# Confirm that exterior border was successfully removed
qtm(US_int_borders)

# Spatial objects, projected into same crs
US <- st_transform(US, crs = Albers_US)
US_int_borders <- US_int_borders %>% st_transform(crs = Albers_US) 

# Function for darkening colors
darken <- function(color, factor=1.11){
  col <- as.vector(col2rgb(color))
  col <- col/factor
  col <- rgb(t(col), maxColorValue=255)
  col
}

# County map, fraction of days with smoke plume coverage --------------------------------------------------------------

# Custom outcome bins and color palette 
bin_cuts <- c(-Inf, 10, 20, 30, 40, Inf)
bin_labels <- c("10", "20", "30", "40", "")

# Utility for blending colors: https://meyerweb.com/eric/tools/color-blend/#:::hex
# Map colors (from light to dark)
palette <- c("#250f55", "#7c288b", "#d54d88", "#f5a8a4", "#faebd3")
# palette <- rev(c("#250f55", "#7c288b", "#d54d88", "#f5a8a4", "#faebd3"))
border_palette <- sapply(palette, darken)
names(border_palette) = NULL

# Map for US 48 

f_smoke_days <- function(name_file, map_year = 0, highlight_county_fips = NA, save_ext = "pdf"){
  # map_year <- 2007
  # highlight_county_fips <- "55025"
  map_year_label <- ifelse(map_year == 0, "2007–2019", map_year)
  
  US_smoked <- US %>% 
    filter(rfrnc_yr == map_year) %>% 
    mutate(smoke_days = cut(hms_deep_1, breaks = bin_cuts)) %>% 
    mutate(highlight = GEOID10 %in% highlight_county_fips)
  all.equal(nrow(US_smoked), 3109) %>% stopifnot()
  table(US_smoked$smoke_days)
  stopifnot(sum(is.na(US_smoked$smoke_days)) == 0)
  
  # Map for US 48 
  m_cont <-
    ggplot(US_smoked) + 
    geom_sf(aes(fill = smoke_days, color = smoke_days), size = 0.05) +
    geom_sf(data = US_int_borders, color = "white", fill = "transparent", size=0.4) +
    geom_sf(data = filter(US_smoked, highlight == 1), color = "white", fill = "transparent", size = 1.2) +
    geom_sf(data = filter(US_smoked, highlight == 1), color = "blue", fill = "transparent", size = 0.6) +
    theme(legend.title=element_blank()) +
    scale_fill_manual(values = rev(palette), drop=FALSE, labels = bin_labels) +    # legend format
    scale_color_manual(values = rev(border_palette), drop=FALSE, labels = bin_labels) +    # legend format
    guides(
      color = "none",
      fill = guide_legend(
        override.aes = list(color = rev(border_palette)), # set color = "transparent" to remove legend key borders
        nrow=1,
        direction = "horizontal",
        title.position = "top",
        label.position = "bottom",
        label.theme = element_text(hjust = 1.3, vjust = 1, size = 11),
        title = paste0("Annual Smoke Days in ", map_year_label))) +
    theme(
      panel.grid.major = element_line(colour = 'transparent'),
      axis.title.x=element_blank(),
      axis.text.x=element_blank(),
      axis.ticks.x=element_blank(),
      axis.title.y=element_blank(),
      axis.text.y=element_blank(),
      axis.ticks.y=element_blank(),
      panel.background=element_blank(),
      panel.grid.minor=element_blank(),
      plot.background=element_blank(),
      legend.title=element_text(face = "bold", size = 12),
      legend.title.align = 0.5,
      legend.position = c(0.72, 0.95),
      legend.key.width = unit(12.5,'mm'),
      legend.key.height = unit(2,'mm'),
      legend.spacing.x = unit(0.6, 'mm'),
      legend.box.background = element_rect(color="white"),
      legend.margin = margin(0,0,0,0),
      plot.margin=unit(c(0,0,0,0),"mm"))
  
  dir.create(file.path(proj_dir, "2_analysis/output_figures"), recursive = TRUE, showWarnings = FALSE)
  if (save_ext == "png") {
    png_res <- 200
    png(file.path(proj_dir, sprintf("2_analysis/output_figures/%s.png", name_file)),
        width = 8.8*png_res, height = 5.3*png_res, pointsize = 4, units = 'px', res = png_res)
  } else if (save_ext == "pdf") {
    pdf(file.path(proj_dir, sprintf("2_analysis/output_figures/%s.pdf", name_file)),
        width =8000/720, height =4900/720)
  }
  
  print(m_cont)
  dev.off()
  
  return(m_cont)
}


# Maps by year
for (year in 2007:2019) {
  f_smoke_days(paste0("figure1_", year), map_year = year, save_ext = "png")
}



# EOF ---------------------------------------------------------------------

